binscatter <- function(x,y,xlab='X',ylab='Y',
                       xlim=c(min(x, na.rm=TRUE),max(x, na.rm=TRUE)),
                       ylim=c(min(y, na.rm=TRUE),max(y, na.rm=TRUE)),
                       pch=19,main = paste0(ylab, ' vs. ', xlab),col = 1,lwd=4){
  
  N <- 20
  
  if (nargs() < 5){
    title = paste0(ylab, ' vs. ', xlab)
  }

  qs <- quantile(x,probs = seq(0,1,1/N))
  ms <- rep(0, length(qs)-1)
  xms <- ms
  for (i in 1:length(qs)-1){
    
    idx <- which(x <= qs[i+1] & x > qs[i])
    ms[i]  <- mean(y[idx])
    xms[i] <- mean(x[idx])
  }
  
  plot(xms,ms, xlab = xlab, ylab = ylab,xlim=xlim,ylim=ylim, main = main, lwd = 3, pch = 19)
  
  mod <- lm(y~x)
  plotx <- c(min(x),max(x))
  plotx <- c(min(xms),max(xms))
  lines(x=plotx,y=predict(mod,newdata=data.frame(x=plotx)),
        col=col,lwd=3)
  
  
#  abline(lm(y ~ x), col = col, lwd = 3)
  
}